function out = RN_OptionPrices(kappa,zeta,eta,sigma,rho_sv,tau,K,r,S,v,X,alpha,PutCall)
                             
%    Returns the option prices

%    kappa      = speed of mean reversion stock return varaince under RN measure
%    zeta       = long-run mean coefficient 0 of stock return varaince under RN measure
%    eta       = long-run mean coefficient 1 of stock return varaince under RN measure
%    sigma      = volatility parameter of stock return variance
%    rho_sv     = correlation between stock return and its variance
    
%    Option Contract features
%    tau        = time to maturity 
%    K          = strike price
%    r          = riskless interest rate
%    S          = current stock price
%    v          = current stock return variance

%    X          = current sentiment
%    alpha      = weight

%   Integration grid
Lphi = 0.000001;  % Lower limit of phi
dphi = 0.01;      % Increment of phi
Uphi = 50;        % Upper limit of phi
phi = [Lphi:dphi:Uphi];
N = length(phi);
% B_OptionPsi(j,phi,kappa_v,v_star_0,sigma_v,rho_sv, tau,K,r,S,v)

%  Computing the conditional probabilities Psi1 and Psi2
int1 = zeros(1,N);
int2 = zeros(1,N);
for k=1:N;
	int1(k) = RN_OptionPsi(kappa,zeta,eta,sigma,rho_sv,tau,K,r,S,v,X,alpha,phi(k),1);
	int2(k) = RN_OptionPsi(kappa,zeta,eta,sigma,rho_sv,tau,K,r,S,v,X,alpha,phi(k),2);
end
I1 = trapz(int1)*dphi;      % trapezoid method for the integral 1 
I2 = trapz(int2)*dphi;      % trapezoid method for the integral 2 
Psi1 = 1/2 + 1/pi*I1;       % conditional probability Psi1 
Psi2 = 1/2 + 1/pi*I2;       % conditional probability Psi2 

CallPrice = S*Psi1 - K*exp(-r*tau)*Psi2;
PutPrice = CallPrice - S + K*exp(-r*tau);

if strcmp(PutCall,'C')
	out = max(CallPrice,0);
else
	out = max(PutPrice,0);
end

